C
C*******************************************************************************
C
C
C                   V A R I A B L E    N A M E    L I S T
C                   _______________    _______    _______
C
C
C  VARIABLE                          DESCRIPTION
C
C  ADT       = TEMPORARY STORAGE VALUE FOR KT/E
C  AII,AJJ   = TEMPORARY STORAGE VALUE FOR DATOS(3,I)
C  AJ        = AREA OF JET NOZZLE (CM**2)
C  ALPA(79)  = KT/E AS INPUT PERMANENT DATA
C  ALPHA(I)  = THERMAL DIFFUSION COEFFICIENT FOR TRACE SPECIES I
C  AM        = MEAN MOLEC. WEIGHT OF GAS MIXTURE(NEGLECTS TRACE SPECIES)
C  APHMIX    = THERMAL DIFFUSIVITY OF MIXTURE AT STATION J
C  APHW      = THERMAL DIFFUSIVITY OF MIXTURE AT STATION W
C  AW        = AREA OF THE COLLECTOR(TARGET), (CM**2)
C
C  C         = CMIX INTERMEDIATE
C  CJ        = DIMENSIONLESS HEAT CAPACITY
C  CMIX      = MIXTURE HEAT CAPACITY AT STATION J (CAL/(G*K))
C  CMH(I)    = RATIO OF EFFECT OF VARIABLE PROP. ON MASS TO HEAT TRANS.
C              OPTIONAL INPUT.LEAVE BLANK IF NO INFORMATION AVAILABLE.
C  COEF1     = COEFFICIENTS FOR COMPUTING HEAT CAPACITY (1000K<T<5000K)
C  COEF2     = COEFFICIENTS FOR COMPUTING HEAT CAPACITY (300K<T<1000K)
C  CONDT     = THERMAL CONDUCTIVITY (TEMPORARY) , (CAL/(CM*SEC*K))
C  CP        = MIXTURE HEAT CAPACITY AT STATION W (CAL/(G*K))
C
C  D         = TEMPORARY STORAGE FOR BINARY DIFF. COEF. (CM**2/SEC)
C  DC        = DISCHARGE COEFFICIENT OF JET NOZZLE. OPTIONAL INPUT.
C              LEAVE BLANK IF NO INFORMATION IS AVAILABLE, IN WHICH
C              CASE THE PROGRAM AUTOMATICALLY ASSIGNS UNITY AS ITS VALUE
C  DA        = DIFF. COEFF. OF TRACE SPECIES AT STATION J (CM**2/SEC)
C  DATOS(1,I)= E/K FOR SPECIES I
C  DATOS(2,I)= SIGMA FOR SPECIES I (COLLISION DIAMETER PARAMETER)
C  DATOS(3,I)= MOLECULAR WEIGHT OF SPECIES I
C  DEN       = DENOMINATOR
C  DIAW      = DIAMETER OF THE CYLINDERICAL COLLECTOR (CM)
C  DIFF      = TEMPORARY STORAGE FOR BINARY DIFF. COEF. (CM**2/SEC)
C  DIV       = DIVERGENCE FACTOR FOR NOZZLE JET DEFINED AS THE RATIO
C              OF CENTERLINE VELOCITY AT COLLECTOR LOCATION TO
C              CENTERLINE VELOCITY AT JET EXIT. DEPENDS ON THE JET
C              DIVERGENCE ANGLE. OPTIONAL INPUT, LEAVE BLANK IF NO
C              INFORMATION AVAILABLE.
C  DJ        = DIAMETER OF JET NOZZLE EXIT PLANE (CM)
C  DW        = DIFF. COEF. OF TRACE SPECIES AT STATION W (CM**2/SEC)
C
C  E/K       = RATIO OF MOLECULAR INTERACTION ENERGY PARAMETER
C              (LENNARD-JONES) TO BOLTZMANN CONSTANT (K)
C  EOKI,EOKJ = TEMPORARY STORAGE FOR DATOS(1,I)=E/K
C  EOVKIJ    = E/K, SQUARE ROOT OF EOKI*EOKJ
C  EPSOVK    = DATOS(1,I)=E/K
C  ERN       = DIFFERENCE BETWEEN PREDICTED AND OBSERVED DEPO. RATES
C  ERROR     = PERCENT ERROR BETWEEN PREDICTED AND OBSERVED DEP. RATE
C  ET        = EXPONENTIAL(-TAU), SORET INTERMEDIATE
C  ETA       = MIXTURE VISCOSITY AT STATION W (G/(CM*SEC))
C  ETAMJ     = TEMPORARY VALUE FOR MIXTURE VISCOSITY (G/(CM*SEC))
C  ETAMIX    = MIXTURE VISCOSITY AT STATION J (G/(CM*SEC))
C
C  F         = FUEL TO AIR RATIO (BY MASS)
C  FSORET    = LOGICAL VALUE PLACED INTO THE INPUT DATA FILE.  IF IT IS
C              SET TO 'T' THEN THE SORET EFFECT WILL BE INCLUDED.  IF IT
C              IS SET TO 'F' THEN THE SORET EFFECT WILL BE EXCLUDED.
C
C  GM        = GAMMA
C  GMR       = (GAMMA-1.)/GAMMA
C  G         = WASSILJEWA-WILKE FACTOR
C
C  ITYPE     = NOT USED BY PROGRAM. OPTIONAL FOR USER PURPOSES.
C  IRUN      = NOT USED BY PROGRAM. OPTIONAL FOR USER PURPOSES.
C
C  LAMDA     = MIXTURE THERMAL CONDUCTIVITY AT STATION W (CAL/(CM*SEC*K))
C  LAMP      = TEMPORARY THERMAL CONDUCTIVITY (CAL/(CM*SEC*K))
C  LAMIX     = MIXTURE THERMAL CONDUCTIVITY AT STATION J (CAL/(CM*SEC*K))
C  LEJ       = LEWIS NUMBER FOR TRACE SPECIES AT STATION J
C  LOD       = RATIO OF MAINSTREAM TURBULENCE LENGTH SCALE TO COLLECTOR DIAM.
C  LW        = LENGTH OF THE CYLINDERICAL COLLECTOR (CM)
C  LEW       = LEWIS NUMBER FOR TRACE SPECIES AT STATION W
C
C  M         = DATOS(3,1)=MOLECULAR WEIGHT OF SPECIES I
C  M(I)      = EQUIVALENT NA2SO4 DEPOSITION RATE OF TRACE SPECIES MASS
C              TRANSFER RATE (MG/HR)
C
C  NU        = NUSSELT NUMBER FOR MASS TRANSFER
C
C  OMEDT,OVDT= COLLISION INTEGRALS
C  OMED(79)  = COLLISION INTEGRAL FOR MASS DIFFUSIVITY
C  OME (79)  = COLLISION INTEGRAL FOR VISCOSITY OR THERMAL CONDUCTIVITY
C
C  PHI       = F(TURB) INTERMEDIATE RELATED TO TURBULENCE INTENSITY
C  PJ        = JET EXIT PLANE PRESSURE (ATM)
C  PO        = STAGNATION PRESSURE (ATM)
C  PR        = PRANDTL NUMBER AT STATION J
C  PROV      = INTERMEDIATE  IN CALCULATION OF ETAMIX
C  PSI       = F(TURB) INTERMEDIATE RELATED TO TURBULENCE LENGTH SCALE
C
C  R         = GAS CONSTANT=1.9872 CAL/(GMOLE*K)
C  RATIO     = SODIUM TO SULFUR MOLAR FLUX RATIO
C  RE        = REYNOLDS NUMBER AT COLLECTOR LOCATION
C  REJ       = REYNOLDS NUMBER AT OUTLET OF JET
C  RHOJ      = DENSITY OF MIXTURE AT JET OUTLET
C  RHOW      = DENSITY OF MIXTURE AT WALL
C
C  SC        = SCHMIDT NUMBER FOR TRACE SPECIES AT STATION J
C  SHAPE     = SHAPE FACTOR (CENTERLINE VELOCITY)/(AVERAGE VELOCITY)
C              OPTIONAL INPUT. LEAVE BLANK IF NO INFORMATION AVAILABLE.
C  SIGMA     = MOLECULAR SIZE PARAMETER (LENNARD-JONES)
C  SIGI,SIGJ = TEMPORARY STORAGE FOR DATOS(2,I)=SIGMA
C  SIGK,SIGL = TEMPORARY STORAGE FOR DATOS(2,I)=SIGMA
C  SIG2      = DATOS(2,I)=SIGMA
C  SIG       = SIGMA VALUE FOR SPECIES I
C  SM(I)     = MASS FLUX OF TRACE SPECIES I (G/CM**2/SEC)
C  SMS       = TOTAL SULFUR MASS FLUX TO THE SURFACE (G/CM**2/SEC)
C  SUM       = TOTAL SODIUM MASS FLUX TO THE SURFACE (G/CM**2/SEC)
C
C  T         = TEMPERATURE (K)
C  TAU       = SORET SUCTION PARAMETER FOR SPECIES I
C  TO        = STAGNATION TEMPERATURE (K)
C  THERM     = THERMAL DIFFUSION FACTOR FOR SPECIES I
C  TJ        = JET (STATIC) TEMPERATURE (K)
C  TURB      = TURBULENCE ENHANCEMENT FACTOR. OPTIONAL INPUT.
C              A) LEAVE BLANK IF NO CONSIDERATION REQUIRED.
C              B) GIVE YOUR OWN ESTIMATE IF INFORMATION AVAILABLE, BUT
C                 LEAVE 'TURIN' AND 'TURL' BLANK IF YOU DO SO.
C  TURIN     = PERCENT MAINSTREAM INTENSITY. OPTIONAL INPUT.
C              A) LEAVE BLANK IF NO CONSIDERATION REQUIRED.
C              B) LEAVE BLANK IF YOU SUPPLY YOUR OWN 'TURB'.
C              C) GIVE YOUR ESTIMATE IF PROGRAM CONSIDERATION DESIRED,
C                 BUT THEN GIVE AN ESTIMATE FOR 'TURL',TOO.
C  TURL      = MAINSTREAM LENGTH SCALE (MACRO) , CM. OPTIONAL INPUT.
C              A) LEAVE BLANK IF NO CONSIDERATION REQUIRED.
C              B) LEAVE BLANK IF YOU SUPPLY YOUR OWN 'TURB'.
C              C) GIVE YOUR ESTIMATE IF PROGRAM CONSIDERATION DESIRED,
C                 BUT THEN GIVE AN ESTIMATE FOR 'TURIN',TOO.
C  TW        = WALL TEMPERATURE (K)
C
C  UJ        = JET EXIT VELOCITY
C
C  VISC(I)   = VISCOSITY OF SPECIES (G/(CM*SEC))
C  V(I)      = NUMBER OF SODIUM ATOMS IN TRACE SPECIES I
C
C  WA        = AIR FLOW RATE (G/SEC)
C  WI,WJ     = TEMPORARY STORAGE VALUES FOR DATOS(3,I)
C  WPRED     = PREDICTED NA2SO4 DEPOSITION RATE (MG/HR)
C  WOBS      = EXPERIMENTALLY OBSERVED NA2SO4 DEPOSTION RATE (MG/HR)
C
C  XJ(I)     = MOLE FRACTION OF SPECIES I IN FLAME
C  X1        = TURIN*RE/100,F(TURB) INTERMEDIATE
C  XW(I)     = MOLE FRACTION OF SPECIES I AT WALL
C  X(4)      = MOLE FRACTION OF DOMINANT SPECIES IN MIXTURE
C
C  YNAJ      = TOTAL SODIUM MOLE FRACTION IN FLAME
C  YNAW      = TOTAL SODIUM MOLE FRACTION AT STATION W
C  YSJ       = TOTAL SULFUR MOLE FRACTION IN FLAME
C  YSW       = TOTAL SULFUR MOLE FRACTION AT STATION W
C
C
C*******************************************************************************
C
C
C  SODIUM VERSION OF CHEMICALLY FROZEN BOUNDARY LAYER THEORY
C  USED FOR SODIUM SULFATE DEPOSITION RATE CALCULATIONS
C
C  MAIN PROGRAM
C
      REAL LAMIX,LAMDA,LEW,LEJ,NU,LW
      REAL M(7)
      DIMENSION V(7),XW(7),XJ(7),SM(7),ALPHA(7),CMH(7)
      COMMON/VA/DATOS(3,5)
      LOGICAL FSORET
      COMMON FSORET
      COMMON/GFR/X(4)
      DATA V/ 1.0 , 1.0 , 2.0 , 0.0 , 0.0 , 1.0 , 0.0 /
C
      WRITE(6,490)
  490 FORMAT(10X,'CFBL THEORY FOR SODIUM SULFATE DEPOSITION RATE',/)
      WRITE(6,900)
  900 FORMAT(10X,'ALL GAS PROPERTIES PERTAIN TO TJ,PJ',/)
C
C**********************************************************************
C                                                                     *
C     PAY SPECIAL ATTENTION TO THE FORMAT OF THE INPUT PARAMETERS     *
C     SEE 'VARIABLE NAME LIST' ABOVE FOR DESCRIPTION OF VARIABLES     *
C                                                                     *
C**********************************************************************
C
C  DOMINANT SPECIES IN THE MIXTURE
C  1=NITROGEN , 2=OXYGEN , 3=WATER , 4=CARBON DIOXIDE
C
C  TRACE SPECIES IN THE MIXTURE
C  1=NAOH , 2=NA , 3=NA2SO4 , 4=SO2 , 5=SO3 , 6=NACL , 7=H2S
C
C
C  INPUT : 1
C
      READ(5,1)IRUN,ITYPE,F,WA,TW,TO,PO,PJ,DC,WOBS,FSORET
    1 FORMAT(2I2,8F9.0,3X,L1)
C
C  INPUT : 2
C
      READ(5,998) DIAW,LW,DJ,TURB,TURIN,TURL,SHAPE,DIV
  998 FORMAT(8F10.0)
C
C  INPUT : 3
C
      DO 502 I=1,7
      READ(5,2) XW(I),XJ(I),CMH(I)
    2 FORMAT(2E15.4,F10.0)
  502 CONTINUE
C
C  OUTPUT OF INPUT PARAMETERS
C
      WRITE(6,978) IRUN,ITYPE,F,TW,TO,PO,PJ,DC,WOBS,FSORET
  978 FORMAT(10X,'IRUN =',I2,/,10X,'ITYPE =',I2,/,10X,'F =',
     *F10.6,/,10X,'TW (K) =',F10.3,/,10X,
     *'TO (K) =',F10.3,/,10X,'PO (ATM) =',F10.6,/,10X,'PJ (ATM) =',F10.6
     *,/,10X,'DC =',F9.6/10X,'WOBS (MG/HR) =',F10.4,/,10X,'FSORET =',L1)
      WRITE(6,988) WA,DIAW,LW,DJ,TURIN,TURL,SHAPE,DIV
  988 FORMAT(10X,'WA (G/SEC) =',F10.5,/,10X,'DIAW (CM) =',F10.5,
     */,10X,'LW (CM) =',F10.5,/,10X,'DJ (CM) =',F10.5,/,10X,
     *'TURIN =',F8.4,'  PERCENT',/,10X,'TURL (CM) =',F10.5,/,10X,
     *'SHAPE =',F7.4,/,10X,'DIV =',F7.4)
      WRITE(6,939)
  939 FORMAT(/,10X,'I=1=NAOH , I=2=NA , I=3=NA2SO4 , I=6=NACL',/)
C
C  ASSIGN E/K,SIGMA AND MOLECULAR WEIGHT OF DOMINANT SPECIES
C
      DATOS(1,1)=71.4
      DATOS(2,1)=3.798
      DATOS(3,1)=28.0134
      DATOS(1,2)=106.7
      DATOS(2,2)=3.467
      DATOS(3,2)=31.9988
      DATOS(1,3)=356.0
      DATOS(2,3)=2.649
      DATOS(3,3)=18.0152
      DATOS(1,4)=195.2
      DATOS(2,4)=3.941
      DATOS(3,4)=44.0098
C
C  THERMAL DIFFUSION FACTOR OF TRACE SPECIES
C  I==NAOH,I=2=NA,I=3=NA2SO4,I=4=SO2,I=5=SO3,I=6=NACL,I=7=H2S
C
      ALPHA(1)=0.1592*(1.-429.5/TW)
      ALPHA(2)=-0.0607*(1.-551.2/TW)
      ALPHA(3)=0.7828*(1.-476.2/TW)
      ALPHA(4)=0.3131*(1.-198.2/TW)
      ALPHA(5)=0.3877*(1.-184.3/TW)
      ALPHA(6)=0.3462*(1.-420.0/TW)
      ALPHA(7)=0.0592*(1.-79.38/TW)
C
C  COMPUTE THE AVERAGE MOLECULAR WEIGHT OF MIXTURE
C  (NEGLECT THE TRACE SPECIES)
C  1=NITROGEN,2=OXYGEN,3=WATER,4=CARBON DIOXIDE
C  MOLE FRACTIONS OF DOMINANT SPECIES AFTER COMBUSTION
C
      DEN=1.+1.0331*F
      X(1)=0.7905/DEN
      X(2)=(0.2095-3.0993*F)/DEN
      X(3)=2.0662*F/DEN
      X(4)=X(3)
      AM=DATOS(3,1)*X(1)+DATOS(3,2)*X(2)+DATOS(3,3)*X(3)+DATOS(3,4)*X(4)
C
C  CALCULATE JET EXIT TEMPERATURE,(K), AND VELOCITY,(CM/SEC)
C
      CALL CALMIX(TO,AM,CMIX)
      GM=CMIX/(CMIX-1.9872/AM)
      GMR=(GM-1.)/GM
      AJ=3.14159*DJ*DJ/4.
      TJ=TO*(PJ/PO)**GMR
      IF(DC.EQ.0.00) DC=1.00
      UJ=((82.057*TJ)/(PJ*AM))*(WA*(1.+F))/(AJ*DC)
C
C  CALCULATE MIXTURE PROPERTIES AT STATIONS J AND W
C
      RHOJ=PJ*AM/TJ/82.057
      RHOW=PJ*AM/TW/82.057
      CALL CALMIX(TJ,AM,CMIX)
      CALL CALMIX(TW,AM,CP)
      CALL MIXPRO(TJ,ETAMIX,LAMIX)
      CALL MIXPRO(TW,ETA,LAMDA)
      APHMIX=LAMIX/CMIX/RHOJ
      APHW=LAMDA/CP/RHOW
      PR=CMIX*ETAMIX/LAMIX
C
C  CALCULATE REYNOLDS NUMBER
C
      IF(SHAPE.EQ.0.0) SHAPE=1.0
      IF(DIV.EQ.0.0) DIV=1.0
      REJ=RHOJ*UJ*DIAW/ETAMIX
      RE=REJ*SHAPE*DIV
C
C  GIVE VALUES OF E/K, SIGMA AND MOLEC. WEIGHT FOR TRACE SPECIES I
C
      DO 10 I=1,7
   25 IF(I.GT.1) GO TO 5
      DATOS(1,5)=1962.0
      DATOS(2,5)=3.804
      DATOS(3,5)=39.9971
      GO TO 11
    5 IF(I.GT.2) GO TO 6
      DATOS(1,5)=1375.
      DATOS(2,5)=3.567
      DATOS(3,5)=22.98977
      GO TO 11
    6 IF(I.GT.3) GO TO 7
      DATOS(1,5)=2221.
      DATOS(2,5)=5.00
      DATOS(3,5)=142.0371
      GO TO 11
    7 IF(I.GT.4) GO TO 8
      DATOS(1,5)=335.4
      DATOS(2,5)=4.112
      DATOS(3,5)=64.0588
      GO TO 11
    8 IF(I.GT.5) GO TO 9
      DATOS(1,5)=431.4
      DATOS(2,5)=4.207
      DATOS(3,5)=80.0582
      GO TO 11
    9 IF(I.GT.6) GO TO 12
      DATOS(1,5)=1989.
      DATOS(2,5)=4.186
      DATOS(3,5)=58.4428
      GO TO 11
   12 CONTINUE
      DATOS(1,5)=301.1
      DATOS(2,5)=3.623
      DATOS(3,5)=34.0758
   11 CONTINUE
C
C  DIFFUSION COEFFICIENT OF TRACE SPECIES AT STATIONS J AND W
C
      CALL DIF(TJ,PJ,DA)
      CALL DIF(TW,PJ,DW)
C
C  LEWIS NUMBER OF TRACE SPECIES AT STATIONS J AND W
C
      LEJ=DA/APHMIX
      LEW=DW/APHW
C
C  SCHMIDT NUMBER OF TRACE SPECIES AT STATION J
C
      SC=ETAMIX/RHOJ/DA
C
C  MASS TRANS. NUSSELT NO.(PERIMETER-AVERAGED,CYLINDERICAL COLLECTOR)
C
      CALL NUM(TW,TO,RE,SC,NU)
C
C  COMPUTE F(SORET)
C
      CALL TMDF(I,TW,TO,LEW,ALPHA,TAU,THERM)
C
C  COMPUTE F(TURB)
C
      CALL TURBL(RE,TURB,TURIN,TURL,DIAW)
C
C  CALCULATE NA2SO4 DEPOSITION RATE CONTRIBUTED FROM SPECIES I
C  BY CFBL FORMULATION.
C
      IF(CMH(I).EQ.0.0) CMH(I)=1.0
      SM(I)=TURB*DA*RHOJ*THERM*((XJ(I)-XW(I))+XW(I)*TAU/THERM*
     *(LEW/LEJ)**0.6*CMIX/CP/CMH(I))*NU*DATOS(3,5)/AM/DIAW
      IF(V(I).EQ.0.0) GO TO 10
      AW=3.14159*DIAW*LW
      M(I)=AW*3.600E6*0.5*142.0371/DATOS(3,5)*SM(I)*V(I)
C
C  PRINT INTERIM OUTPUT
C
      IF(I.GT.1) GO TO 546
      IF(FSORET) GO TO 100
      WRITE(6,544)
  544 FORMAT(/,10X,'CALCULATIONS DO NOT INCLUDE THERMAL DIFFUSION',/)
  100 WRITE(6,545)
  545 FORMAT(T2,'I',T11,'D(I)',T25,'NU(I)',T40,'SC(I)',T55,'X(I),
     CJ',T70,'X(I),W',T85,'TAU(I)',T100,'F(SORET),I',T116,'M(I)')
  546 WRITE(6,550) I,DA,NU,SC,XJ(I),XW(I),TAU,THERM,M(I)
  550 FORMAT(1X,I1,6E15.5,F15.5,2E15.5,/)
   10 CONTINUE
C
C  COMPUTE TOTAL SODIUM(NA) AND SULFUR(S) MASS FLUX
C     SUM=J(NA), SMS=J(S)
C
      SUM=SM(1)*22.98977/39.9971+SM(2)*1.0+SM(3)*2.*22.98977
     //142.0371+SM(6)*22.98977/58.4428
      SMS=SM(3)*32.06/142.0371+SM(4)*32.06/64.0588+SM(5)*
     *32.06/80.0582+SM(7)*32.06/34.0758
C
C  COMPUTE SODIUM TO SULFUR MOLAR FLUX RATIO
C
      RATIO=SUM/SMS*32.06/22.98977
C
C  COMPUTE TOTAL NA AND S MOLE FRACTIONS AT STATIONS J AND W
C
      YNAJ=XJ(1)+XJ(2)+(XJ(3)*2)+XJ(6)
      YNAW=XW(1)+XW(2)+(XW(3)*2)+XW(6)
      YSJ=XJ(3)+XJ(4)+XJ(5)+XJ(7)
      YSW=XW(3)+XW(4)+XW(5)+XW(7)
      WRITE(6,211)
  211 FORMAT(//,T8,'YNAJ',T24,'YNAW',T40,'YSJ',T57,'YSW')
      WRITE(6,212)YNAJ,YNAW,YSJ,YSW
  212 FORMAT(4E16.6,/)
C
C  CALCULATE NA2SO4 DEPOSITION RATE (MG/HR)
C
      WPRED=M(1)+M(2)+M(3)+M(6)
C
C  COMPUTE PERCENT ERROR BETWEEN M(PRED) AND M(OBS)
C
      ERR=WPRED-WOBS
      ERROR=ERR/WOBS*100.
C
C  PRINT RESULTS
C
      WRITE(6,4)AM,TURB,RE,SUM,SMS,RATIO,WPRED,ERROR,GM,PR,X(1),X(2),
     *X(3),X(4),TJ
    4 FORMAT(/,10X,'AM=',F10.5,/,10X,'TURB=',F10.5,/,10X,'REEFF=',F10.3,/,
     *10X,'SUM (G/CM**2/SEC) =',E13.5,/,10X,'SMS (G/CM**2/SEC) =',E13.5,/,
     *10X,'RATIO =',F10.4,/,10X,
     *'PREDICTED NA2SO4 DEPOSITION RATE (MG/HR) =',F10.5,/,10X,
     *'ERROR (%) =',F10.4,/,10X,'GM =',F10.5,/,10X,'PR =',F10.6,/,10X,
     *'X(N2) =',F10.6,/,10X,'X(O2) =',F10.6,/,10X,'X(H20) =',F10.6,/,10X,
     *'X(CO2) =',F10.6,/,10X,'TJ (K) =',F10.4,/)
      WRITE(6,950)
  950 FORMAT(/,T5,'RHOJ(G/CM**3)',T21,'UJ(CM/S)',T37,'ETAMIX(POISE)
     C',T52,'LAMIX(CAL/CM/K/S)',T73,'CMIX(CAL/G/S)')
      WRITE(6,951)RHOJ,UJ,ETAMIX,LAMIX,CMIX
  951 FORMAT(4E16.6,4X,E16.6,//)
      STOP
      END
C
C*******************************************************************************
C
      SUBROUTINE PARAM (J,EPSOVK,SIG2,M)
C
C  SUBROUTINE PARAM GIVES THE VALUES OF E/K,SIGMA AND MOLECULAR WEIGHT
C
      REAL M
      COMMON/VA/DATOS(3,5)
      EPSOVK=DATOS(1,J)
      SIG2=DATOS(2,J)
      M=DATOS(3,J)
      RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE CALESL(J,T,CJ)
C
C  THIS SUBROUTINE CALCULATES DIMENSIONLESS HEAT CAPACITY FROM A
C  CURVE FIT VALID FOR  300K < T < 1000K AND 1000K < T < 5000K
C
      DIMENSION COEF1(5,4),COEF2(5,4)
      DATA COEF1/ 0.28532899E+01 , 0.16022128E-02 ,-0.62936893E-06 ,
     *            0.11441022E-09 ,-0.78057465E-14 ,
     *            0.36122139E+01 , 0.74853166E-03 ,-0.19820647E-06 ,
     *            0.33749008E-10 ,-0.23907374E-14 ,
     *            0.26340654E+01 , 0.31121899E-02 ,-0.90278449E-06 ,
     *            0.12673054E-09 ,-0.69164732E-14 ,
     *            0.44608041E+01 , 0.30981719E-02 ,-0.12392571E-05 ,
     *            0.22741325E-09 ,-0.15525954E-13 /
      DATA COEF2/ 0.37044177E+01 ,-0.14218753E-02 , 0.28670392E-05 ,
     *           -0.12028885E-08 ,-0.13954677E-13 ,
     *            0.37837135E+01 ,-0.30233634E-02 , 0.99492751E-05 ,
     *           -0.98189101E-08 , 0.33031825E-11 ,
     *            0.41675564E+01 ,-0.18106868E-02 , 0.59450878E-05 ,
     *           -0.48670871E-08 , 0.15284144E-11 ,
     *            0.24007797E+01 , 0.87350957E-02 ,-0.66070878E-05 ,
     *            0.20021861E-08 , 0.63274039E-15 /
      IF(T.GT.1000.) GO TO 10
      CJ=COEF2(5,J)*T**4+COEF2(4,J)*T**3+COEF2(3,J)*T**2+COEF2(2,J)*T
     *+COEF2(1,J)
      GO TO 20
   10 CJ=COEF1(5,J)*T**4+COEF1(4,J)*T**3+COEF1(3,J)*T**2+COEF1(2,J)*T
     *+COEF1(1,J)
   20 RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE CALMIX(T,AM,CMIX)
C
C SUBROUTINE CALMIX PROVIDES MIXTURE HEAT CAPACITY IN CAL/(G*K)
C
      COMMON/GFR/X(4)
      C=0.
      DO 8 J=1,4
      CALL CALESL(J,T,CJ)
      C=C+CJ*X(J)*1.9872
    8 CONTINUE
      CMIX=C/AM
      RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE COLINT(II,JJ,T,OVDT,ETA,CONDT,OMEDT)
C
C  THIS SUBROUTINE GIVES COLLISION INTEGRALS,VISCOSITY AND THERMAL
C  CONDUCTIVITY OF SPECIES USING CHAPMAN-ENSKOG THEORY
C
      REAL M
      DIMENSION OME(79),ALPA(79),OMED(79)
      DATA ALPA/ 0.300 , 0.350 , 0.400 , 0.450 , 0.500 , 0.550 , 0.600 ,
     *           0.650 , 0.700 , 0.750 , 0.800 , 0.850 , 0.900 , 0.950 ,
     *           1.000 , 1.050 , 1.100 , 1.150 , 1.200 , 1.250 , 1.300 ,
     *           1.350 , 1.400 , 1.450 , 1.500 , 1.550 , 1.600 , 1.650 ,
     *           1.700 , 1.750 , 1.800 , 1.850 , 1.900 , 1.950 , 2.000 ,
     *           2.100 , 2.200 , 2.300 , 2.400 , 2.500 , 2.600 , 2.700 ,
     *           2.800 , 2.900 , 3.000 , 3.100 , 3.200 , 3.300 , 3.400 ,
     *           3.500 , 3.600 , 3.700 , 3.800 , 3.900 , 4.000 , 4.100 ,
     *           4.200 , 4.300 , 4.400 , 4.500 , 4.600 , 4.700 , 4.800 ,
     *           4.900 , 5.000 , 6.000 , 7.000 , 8.000 , 9.000 , 10.00 ,
     *           20.00 , 30.00 , 40.00 , 50.00 , 60.00 , 70.00 , 80.00 ,
     *           90.00 , 100.0 /
      DATA OME/ 2.7850 , 2.6280 , 2.4920 , 2.3680 , 2.2570 , 2.1560 ,
     *          2.0650 , 1.9820 , 1.9080 , 1.8410 , 1.7800 , 1.7250 ,
     *          1.6750 , 1.6290 , 1.5870 , 1.5490 , 1.5140 , 1.4820 ,
     *          1.4520 , 1.4240 , 1.3990 , 1.3750 , 1.3530 , 1.3330 ,
     *          1.3140 , 1.2960 , 1.2790 , 1.2640 , 1.2480 , 1.2340 ,
     *          1.2210 , 1.2090 , 1.1970 , 1.1860 , 1.1750 , 1.1560 ,
     *          1.1380 , 1.1220 , 1.1070 , 1.0930 , 1.0810 , 1.0690 ,
     *          1.0580 , 1.0480 , 1.0390 , 1.0300 , 1.0220 , 1.0140 ,
     *          1.0070 , 0.9999 , 0.9932 , 0.9870 , 0.9811 , 0.9755 ,
     *          0.9700 , 0.9649 , 0.9600 , 0.9553 , 0.9507 , 0.9464 ,
     *          0.9422 , 0.9382 , 0.9343 , 0.9305 , 0.9269 , 0.8963 ,
     *          0.8727 , 0.8538 , 0.8379 , 0.8242 , 0.7432 , 0.7005 ,
     *          0.6718 , 0.6504 , 0.6335 , 0.6194 , 0.6076 , 0.5973 ,
     *          0.5882 /
      DATA OMED/ 2.6620 , 2.4760 , 2.3180 , 2.1840 , 2.0660 , 1.9660 ,
     *           1.8770 , 1.7980 , 1.7290 , 1.6670 , 1.6120 , 1.5620 ,
     *           1.5170 , 1.4760 , 1.4390 , 1.4060 , 1.3750 , 1.3460 ,
     *           1.3200 , 1.2960 , 1.2730 , 1.2530 , 1.2330 , 1.2150 ,
     *           1.1980 , 1.1820 , 1.1670 , 1.1530 , 1.1400 , 1.1280 ,
     *           1.1160 , 1.1050 , 1.0940 , 1.0840 , 1.0750 , 1.0570 ,
     *           1.0410 , 1.0260 , 1.0120 , 0.9996 , 0.9878 , 0.9770 ,
     *           0.9672 , 0.9576 , 0.9490 , 0.9406 , 0.9328 , 0.9256 ,
     *           0.9186 , 0.9120 , 0.9058 , 0.8998 , 0.8942 , 0.8888 ,
     *           0.8836 , 0.8788 , 0.8740 , 0.8694 , 0.8652 , 0.8610 ,
     *           0.8568 , 0.8530 , 0.8492 , 0.8456 , 0.8422 , 0.8124 ,
     *           0.7896 , 0.7712 , 0.7556 , 0.7424 , 0.6640 , 0.6232 ,
     *           0.5960 , 0.5756 , 0.5596 , 0.5464 , 0.5352 , 0.5256 ,
     *           0.5170 /
      IF(JJ.EQ.0)GO TO 67
      CALL PARAM(JJ,EOKJ,SIGL,AJJ)
      CALL PARAM(II,EOKI,SIGK,AII)
      M=(AII+AJJ)/2.
      SIG2=(SIGL+SIGK)/2.
      EOVKIJ=SQRT(EOKJ*EOKI)
      ADT=T/EOVKIJ
      GO TO 68
   67 CALL PARAM(II,EPSOVK,SIG2,M)
      ADT=T/EPSOVK
   68 J=3
   69 IF(ADT-ALPA(J))112,111,110
  110 J=J+1
      IF(J.LE.79)GO TO 69
      GO TO 204
  111 OVDT=OME(J)
      OMEDT=OMED(J)
      GO TO 200
  112 B=(ADT-ALPA(J-1))/(ALPA(J)-ALPA(J-1))
      OVDT=B*(B-1.)*OME(J-2)/2.+(1.-B**2)*OME(J-1)+B*(B+1.)*OME(J)/2.
      OMEDT=B*(B-1.)*OMED(J-2)/2.+(1.-B**2)*OMED(J-1)+B*(B+1.)*OMED(J)/2
     *.
      GO TO 200
  204 OVDT=0.5882*((ADT/100.)**(-0.17))
      OMEDT=0.5170*((ADT/100.)**(-0.15))
  200 ETA=26.693*SQRT(M*T)/((SIG2)**2*OVDT)*1.E-06
      CALL CALESL(II,T,CR)
      CONDT=(1.9872/M)*(15./4.+1.32*(CR-5./2.))*ETA
      RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE MIXPRO(T,ETAMIX,LAMIX)
C
C  THIS SUBROUTINE PROVIDES MIXTURE VISCOSITY AND THERMAL CONDUCTIVITY
C
      REAL LAMP,LAMIX
      DIMENSION VISC(4),LAMP(4)
      COMMON/VA/DATOS(3,5)
      COMMON/GFR/X(4)
      ETAMIX=0.
      LAMIX=0.
      I1=0
      DO 40 J=1,4
      N=J
      CALL COLINT(N,I1,T,OVDT,ETAMJ,CONDT,OMEDT)
      VISC(J)=ETAMJ
      LAMP(J)=CONDT
   40 CONTINUE
      DO 18 I=1,4
      PROV=0.
      DO 28 J=1,4
      G=1./(2.*SQRT(2.))*((1.+DATOS(3,I)/DATOS(3,J))**(-0.5))*(1.
     * +SQRT(VISC(I)/VISC(J))*SQRT(SQRT(DATOS(3,J)/DATOS(3,I))))**2
      PROV=PROV+G*X(J)
   28 CONTINUE
      ETAMIX=ETAMIX+X(I)*VISC(I)/PROV
      LAMIX=LAMIX+X(I)*LAMP(I)/PROV
   18 CONTINUE
      RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE COEDIF(T,P,I,J,D)
C
C  THIS SUBROUTINE CALCULATES BINARY DIFFUSION COEFFICIENT
C
      CALL COLINT(I,J,T,OVDT,ETA,CONDT,OMEDT)
      CALL PARAM(I,EPSOVK,SIGI,WI)
      CALL PARAM(J,EPSOVK,SIGJ,WJ)
      SIG=1./2.*(SIGI+SIGJ)
      D=0.0018583*SQRT((T**3)*(1./WI+1./WJ))/(P*OMEDT*(SIG**2))
      RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE DIF(T,P,DA)
C
C  SUBROUTINE DIF GIVES DIFFUSION COEFFICIENT OF TRACE SPECIES K
C  IN GAS MIXTURE
C
      COMMON/GFR/X(4)
      DEN=0.0
      DO 101 J=1,4
      K=5
      CALL COEDIF(T,P,K,J,DIFF)
      DEN=DEN+X(J)/DIFF
  101 CONTINUE
      DA=1./DEN
      RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE TMDF(I,TW,TO,LEW,ALPHA,TAU,THERM)
C
C  SUBROUTINE TMDF CALCULATES THERMAL DIFFUSION EFFECT, F(SORET),
C  LABELLED AS 'THERM'.  IF FSORET IS FALSE, THERM=1.0 .
C
      LOGICAL FSORET
      COMMON FSORET
      REAL LEW
      DIMENSION ALPHA(7)
      IF(FSORET) GO TO 10
      THERM=1.0
      TAU=0.0
      GO TO 20
   10 TAU=ALPHA(I)*LEW**0.4*(TO-TW)/TW
      ET=EXP(-TAU)
      THERM=TAU/(1.-ET)
   20 RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE NUM(TW,TO,RE,SC,NU)
C
C  SUBROUTINE NUM CALCULATES PERIMETER AVERAGED NUSSELT NUMBER
C  INCLUDING EFFECTS OF MACH NUMBER AND VARIABLE PROPERTIES
C
      REAL NU
      NU=(0.40*SQRT(RE)+0.06*RE**(2./3.))*(SC**0.4)*
     1(TO/TW)**(-0.04)
      RETURN
      END
C
C*******************************************************************************
C
      SUBROUTINE TURBL(RE,TURB,TURIN,TURL,DIAW)
C
C  SUBROUTINE TURBL CALCULATES EFFECT OF MAINSTREAM TURBULENCE
C
      REAL LOD
      IF(TURIN.NE.0.0) GO TO 5
      IF(TURB.EQ.0.0) GO TO 10
      GO TO 30
    5 LOD=TURL/DIAW
      IF(LOD.GT.2.0) PSI=0.124*(LOD-11.0)**2.+2.
      IF(LOD.LE.2.0) PSI=-4.0*(LOD-1.75)**2.+12.25
      PSI=PSI/1000.
      XI=TURIN*RE/100.
      IF(XI-100.)1,2,2
    1 PHI=12.375*(1.-(1.-XI/100.)**1.5)
      GO TO 20
    2 PHI=9.0+0.03375*XI
   20 TURB=1.+PSI*PHI
      GO TO 30
   10 TURB=1.0
   30 RETURN
      END
